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ABSTRACT 

We consider two phase accretion disk-corona models for active galactic nuclei and 
some X-ray binaries. We describe in detail how one can exactly solve the polarized 
radiative transfer and Comptonization using the iterative scattering method, while 
simultaneously solving the energy and pair balance equation for both the cold and hot 
phases. We take into account Compton scattering, photon-photon pair production, 
pair annihilation, bremsstrahlung, and double Compton scattering, as well as exact 
reflection from the cold disk. We consider coronae having slab geometry as well 
as coronae consisting of one or more well separated active regions of cylinder or 
hemisphere geometry. 

The method is useful for determining the spectral intensity and the polarization 
emerging in different directions from disk-corona systems. The code is tested against 
a Monte-Carlo code. We also compare with earlier, less accurate, work. The method 
is more than an order of magnitude faster than applying Monte Carlo methods to the 
same problem and has the potential of being used in spectral fitting software such as 
XSPEC. 



Subject headings: accretion disks - galaxies: Seyfert - gamma rays: theory 
polarization - radiative transfer - X-rays: general 



1. Introduction 

Both active galactic nuclei (AGN) and certain X-ray binaries (the galactic black hole 
candidates, GBHC) show X-ray spectra extending into the hard X-rays (e.g. Mushotzky, Done, & 
Pounds 1993; Tanaka & Levin 1995). The X-ray spectra of Seyfert 1 galaxies show at least two 
components: 1) an intrinsic power law component with an intensity index, a ~ 0.9 — 1.0, in the 
2-18 ke V range and a spectral cutoff at a few hundred keV (Zdziarski et al. 1994; Madejski et al. 
1995; Zdziarski et al. 1995), and 2) a superimposed reflection component arising from reflection 
and reprocessing of the intrinsic power law by cold opaque matter subtending ~ 1-2 it solid angle 
as viewed from the X-ray source (e.g. Nandra & Pounds 1994). GBHC such as Cyg X-l and 
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1E1740. 7-2942 show power law spectra extending up to hundreds of keV (e.g. Gilfanov et al. 1994; 
Tanaka & Levin 1995). The characteristic features of reflection have been seen in GBHC as well 
(e.g. Done et al. 1992). 

AGN and GBHC are believed to be powered by accretion through an accretion disk. In the 
unified model for AGN (e.g. Antonucci 1993), it is believed that we are viewing the disks in 
Seyfert 1 galaxies more or less face on, while for Seyfert 2 galaxies we are viewing the disk more or 
less edge on through a molecular torus. In GBHC sources we are viewing the binary system along 
some given (possibly time dependent) direction. The X-ray spectra indicate the existence of both 
hot X-ray emitting and cold reflecting gas components. The exact geometry is not known, but a 
currently popular model is the two-phase disk-corona model (e.g. Haardt & Maraschi 1991, 1993, 
hereafter HM93). The black body disk radiation from the cold disk (in the EUV for AGN, and in 
the soft X-rays for GBHC) enters the hot corona from one side and gets Comptonized into the 
X-rays. Part of this X-ray radiation is incident on the cold disk and is partly reflected but mainly 
reprocessed into soft black body radiation. The remaining part forms the X-ray spectrum leaving 
the disk-corona system. Both the black body and the Comptonized spectrum are anisotropic 
so observers at different viewing angles see different spectra (HM93). It immediately clear that 
such models are neither homogeneous nor spherically symmetric. In order to correctly interpret 
observed X-ray spectra of AGN and GBHC one needs to know the theoretical spectra for mildly 
relativistic temperature and for different viewing angles. 

Theoretical Comptonized spectra have been computed for two decades now. Almost all 
work make simplifying assumptions that render them useless for interpreting X-ray spectra from 
sources were anisotropic effects are important. We briefly discuss standard methods for modeling 
Comptonized spectra from AGN and GBHC at mildly relativistic temperatures. 

One approach is to treat the photon and pair producing processes as well as the energy and 
pair balance in great detail, but to make large simplifications regarding the radiative transfer 
using various prescriptions for the spectral shape and using simple escape probabilities to get the 
photon density (e.g. Zdziarski 1985; Pietrini & Krolik 1995). Such calculations can only give very 
approximate relations between the typical spectral shape and other parameters. 

Other approaches are to do detailed radiative transfer using Monte Carlo methods for 
geometries such as slabs or spheres (e.g. Hua & Titarchuk 1995) or to improve the analytical 
theory of Comptonization (e.g. Titarchuk 1994). Normally the following simplifications are made: 
1) processes other than Comptonization are neglected, 2) pair balance is not imposed, 3) angle 
dependence of output spectra is not considered, 4) the soft photon injection is homogeneous 
throughout the source. There are some exceptions. For example, Skibo et al. (1995) includes 
bremsstrahlung, pair production, and pair balance, and Zdziarski et al. (1994) assume the soft 
photons to be injected at one of the slab surfaces. Most work neglect reflection by cold matter. 
The few papers considering polarized radiative transfer (e.g. Sunyaev & Titarchuk 1985; Haardt 
& Matt 1993) use the Rayleigh matrix, which is not valid for temperatures and photon energies 
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above 50 keV. 

Pioneering steps were taken by Haardt & Maraschi (1991, 1993) and Haardt (1993, 1994) to 
solve for angle-dependent spectra from disk-corona systems including reflection as well as energy 
and pair balance. The method treated the first scattering order accurately, but assumed higher 
order scatterings to be isotropic. Compton recoil was neglected and therefore the spectral cutoff 
at photon energies around and above kT e could not be treated. Furthermore, the pair balance 
treatment was approximate as they adopted the semi-analytical theory of Zdziarski (1985) using 
prescribed spectra. Finally, when determining the reflected spectrum they assumed the X-ray 
intensity incident on the cold matter to be isotropic. Further steps were taken by Stern et al. 
(1995a, b) who used nonlinear Monte Carlo techniques to treat disk-corona systems with different 
inhomogeneous coronal geometries. Here, exact pair balance was imposed throughout the coronal 
region. Furthermore, output spectra (including the cutoff) had sufficiently good statistics to allow 
the dependence on the viewing angle to be studied. These work found that anisotropic effects are 
very important for spectra from disk-corona systems. At mildly relativistic temperatures the first 
order scattering of the soft disk radiation is suppressed in the face-on direction. This causes the 
face-on spectrum to be harder up to the spectral peak of the second order scattering, where an 
anisotropy break appears. Above the break the Comptonized spectrum in all directions resembles 
the angle-averaged one. As shown in Stern et al. (1995b), the anisotropy break easily appears in 
the 2-18 keV range making a.i— xs strongly dependent on viewing angle. The reflection component 
is also strongly anisotropic (e.g. Matt 1993; Magdziarz & Zdziarski 1995; Poutanen, Nagendra, & 
Svensson 1996). 

The disadvantages with the nonlinear Monte Carlo method is that it is computer intensive 
requiring more than one hour per run on a Sparc 20. In order to be able to do spectral fitting of 
observed spectra one needs a much faster method to compute accurate angle dependent spectra 
from disk-corona systems. In the present paper we describe such a fast code based on the iterative 
scattering method, i.e. where the radiative transfer equation is solved for each scattering order 
(e.g. Sunyaev & Titarchuk 1985). Most important radiation processes and photon-photon pair 
production are included. Energy and pair balance can be imposed. The reflection is accurately 
treated accounting for the full angular dependence of the incident spectrum (Poutanen et al. 
1996). Both slab, cylinder, and hemisphere geometries of the corona can be treated. Finally, 
the radiative transfer is polarized, both as regards the Comptonization and the reflection. The 
code has already been used to interpret the statistics of observed X-ray spectral indices and 
compactnesses from Seyfert 1 galaxies (Stern et al. 1995b). The purpose of the present paper is 
to fully document the methods used in the code in a selfcontained way. 

In the remainder of the paper, we first describe the setup of the two-phase disk-corona model 
in § The methods of solving the radiative transfer equation in different geometries are considered 
in § ||. The energy and pair balance and details of the iteration procedure are considered in § 
We compare our results with those of other available codes in § ||, where we also consider the 
accuracy and efficiency of various approximations that can be employed in order to decrease the 
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computing time. Finally, we summarize our work in § ||. Expressions for the reaction rates and 
redistribution functions (i.e. the Compton redistribution matrix) are given in the Appendices. 

2. Setup 

We consider the simple two-phase disk-corona model, where a hot corona is located above 
an optically thick plane-parallel cold slab ("disk") (e.g. Haardt & Maraschi 1991, 1993; Haardt, 
Maraschi, & Ghisellini 1994). The hot corona is either a plane-parallel slab (with vertical thickness 
H), or an active region with the shape of a hemisphere (with radius R), or a cylinder (with 
vertical height, H, and horizontal radius, R). We allow energy dissipation in both the corona 
and the cold disk. The radiation escaping from the cold disk consists of a soft component, and 
a reflected component. The soft flux is equal to the sum of the absorbed incident flux from the 
corona, and the flux due to local energy dissipation in the cold disk. The spectral shape of these 
soft components is assumed to be Planckian with temperatures Tbb and Tdi s k> respectively (note 
that Tbb > ?disk)- The shape of the reflected component is determined by the shape of incident 
coronal X-ray radiation (mainly resulting from Comptonization of the soft disk radiation) and 
the effects of photoelectric absorption and Compton scattering in the cold disk (see, e.g., White, 
Lightman, & Zdziarski 1988; Magdziarz & Zdziarski 1995). 

We can treat coronae both with or without pairs. In this paper, we mostly consider the case 
of a pure pair corona without any background plasma. The pure pair corona is a consequence of 
photon-photon pair production above the cold disk. The electrons and positrons are assumed to 
have a relativistic Maxwellian distribution of the same temperature, = £T e /m e c 2 . The corona is 
assumed to be uniform in temperature and pair density, and pair escape is neglected. If all power is 
dissipated in the corona then, for a given geometry, the two parameters: the total power dissipated 
in corona, Ldiss> and the temperature, Tbb, of the reprocessed radiation, uniquely determine the 
optical depth, tt, and the coronal temperature, 0. If the radiation produced internally in the disk 
is not negligible, then two more parameters are important: the disk temperature, Tjisk, and the 
ratio, d = ^disk/Aiiss) where Ldisk is the luminosity that is produced internally in the disk and 
that enters the corona. For all geometries, tt is defined as the total vertical Thomson optical 
depth of the corona (along the symmetry axis in the case of hemisphere geometry). 

To solve the pair balance equation, the energy balance equations for the cold and hot phases, 
and the radiative transfer in the corona self-consistently, we make use of an iteration procedure. 
To reduce computing time we choose to fix 0, which allow us to compute the thermal Compton 
redistribution matrix and cross section, and the coronal emissivities for pair annihilation and 
bremsstrahlung before doing the iterations. We then adjust tt and Tdi ss until the radiation spectra 
from solving the radiative transfer satisfy the energy balance equations and the pair balance. 

When solving the radiative transfer/Comptonization problem, we account for the angular 
anisotropy of the radiation as well as its polarization properties. The reprocessing in the cold 
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disk is described by Green's matrix (consisting of four Green's functions) for reflection, where we 
fully account for the Compton effect, photoelectric absorption, and iron fluorescence, as well as 
for the angular and polarization properties of the radiation incident on the cold disk (Poutanen 
et al. 1996). We reduce the time needed to compute the reflection spectrum substantially using 
the precalculated Green's matrix and achieve better accuracy than all previous treatments of the 
problem (e.g., White et al. 1988; Magdziarz & Zdziarski 1995). 

As additional photon sources and cooling processes of the corona, we consider electron- 
electron, positron-positron, and electron-positron bremsstrahlung, double Compton scattering, 
and pair annihilation. We also account for photon absorption due to pair production, which can 
be important in determining the spectral shapes at pair producing energies (hu > m e c 2 ) and, thus, 
in influencing the pair balance. 

The radiative transfer equation is solved by expanding the radiation field in scattering orders 
(the iterative scattering method, e.g. Sunyaev & Titarchuk 1985). The intensity in a slab-corona 
is a function of vertical position (i.e. the Thomson optical depth variable, r), zenith angle and 
frequency but is azimuth- independent. In an active region, the intensity, of course, depends on the 
distance from the symmetry axis and the azimuth angle making the spatial part of the problem 
two-dimensional. However, by averaging the radiation field over horizontal layers in the active 
region, we convert the 2D-problem into a ID-problem suitable for our ID-code. We discuss the 
accuracy of this conversion in § |j| 

In the pair balance, we use a volume-averaged pair production rate. In the energy balance 
equations, we need the total luminosities emerging from the cold and the hot phases. Therefore 
we compute and sum the radiative fluxes emerging from all surfaces of the disk and the corona 
accounting for all radiative transfer effects. 



3. Radiative Transfer 

3.1. Radiative Transfer in a Slab Corona 

Due to azimuthal symmetry and the absence of sources of circular polarization, the radiation 
field and the degree of polarization at vertical position z can be fully described by a Stokes vector 
consisting of two Stokes parameters (Chandrasekhar 1960) / = I(z,x,fi) = (I, Q) T , where T 
denotes the transposed vector. The radiative transfer equation describing the propagation of 
polarized light through a plane-parallel electron (and positron) atmosphere in steady-state can be 
written in the following form: 

dl(z x,fi) = -(n e(7cB (x) + a-fJz, x, n))I{z, x, //) + n e (PrS(z, x, fi) + e(x), (1) 
dz 

where x = hv/rrieC 2 is the photon energy, n e = n + +n_ is the total electron and positron density, fj, is 
the cosine angle between the slab normal and the direction of photon propagation, cr cs (x) cm 2 is the 
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thermal Compton scattering cross section, cjt cm 2 is the Thomson cross section, a 77 (z,x, / u) cm -1 
is the absorption coefficient due to photon-photon pair production, and S(z, x, /x) is the electron 

scattering source function. The emissivity, e(x) = e ++ + e + e H + e ann + enc erg cm -1 s _1 sr _1 , 

which is assumed to be isotropic and homogeneous, includes all photon sources in the atmosphere, 
in our case electron-electron, positron-positron, electron-positron bremsstrahlung, annihilation, 
and double Compton radiation. The expressions for the emissivities and absorption coefficients 
are given in the Appendices. 

Using the following notation for the dimensionless intensity, source function, emissivities, 
scattering cross section and absorption coefficient: 



I' = I ^ S' = S 



e = e r, <7~, = — , a. 



77 



Tn e c m e c? n e m e cy or n e aT 

and removing the primes, the radiative transfer equation can be written in the following 
dimensionless form: 

^ d/(r,x,^) = -(a CB (x) + a 77 (r, x, (i))I(t, x, /i) + S{t, x, /i) + e(x), (2) 

where dr = ain e dz is the differential Thomson optical depth. Hereafter we will only use 
dimensionless quantities (except in the Appendices). 

The Thomson optical depth of the slab is tt = H<j-xn e . The boundary conditions at the 
upper and lower surface of the slab are: 

I{t = 7"Ti x, — n) = 0, \i > 0, 

l(r = 0,x,/x) = I in (x,fi), (3) 

i.e. there is no radiation incident at the upper surface, and the radiation incident at the lower 
surface consists of a reflected component, a soft reprocessed component, and a soft component 
internally produced in the disk: 

L m (x, fl) = I rcfl (x, H) + C hh B x (T hh ) + c disk B x (T Aisk ). (4) 

The reflected radiation, I Te ^(x, /j), from the cold disk can be found as a convolution of a reflection 
matrix (Green's matrix) G(x, //; x%, fj,x) with the incident radiation: 

_ roc rl _ 

I Te& (x,n)= dxi / dfj,iG(x, fj,; x 1} Hx)I{r = Q,xx,-fJ,x)- (5) 
Jx Jo 

The Green's matrix maps incoming radiation at (xx, —fJ-x) into reflected radiation at (x, /i). 
To compute Green's matrix we use the method developed by Poutanen et al. (1996). The 
normalization constants in front of the Planckian functions in equation (|j) are determined by 
normalizing the black-body flux to the soft compactnesses, l vepr and l d [ sk (to be defined in § |Q| ), 
as: 

POO POO 

c hh ir / dxB x (T hh ) = Z ropr ; Cdisk^r / dxB x (T disk ) = l disk . (6) 
Jo Jo 
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The thermal electron scattering source function, S(t,x,(j,), can be expressed in terms of the 
azimuth-averaged Compton redistribution matrix (see e.g. Poutanen &; Vilhu 1993): 

S(t,x,h)=x / — o- / d//i J(r,xi,^i), (7) 



X| 7_i ^ fi 2 l i?22 

or in operator form as 

5 = RJ . (8) 

The factor x 1 jx\ in equation @ appears because we use the photon intensity instead of the photon 
occupation number to describe the radiation field (see Nagirner & Poutanen 1994). Expressions 



for the redistribution functions R are given in Appendix A.l 



We solve the integro-differential equation (||) by expanding the Stokes vector / in scattering 
orders (Neumann series): 



* = £**, (9) 

k=0 

where 1^ is the Stokes vector for photons having undergone k scatterings (see, e.g., Sunyaev & 
Titarchuk 1985). This expansion converges quickly for sufficiently small optical depths (tt Ss 1). 
The source function for the non-scattered component consist of the flux incident on the corona at 
the bottom surface, r = 0, and of the internal coronal sources: 

S (t, x, (jt) = fj,l in (x, h)H(h)5(t) + e(x) (10) 

where H(y) is the Heaviside function. The Stokes vectors, I^(r,x,fi) = Ifc(r, x, ±/x), for the 
upward and downward radiation and for all scattering orders k > are calculated employing the 
iteration formulae: 

r dr' ~ { r dr" 1 

IZ(T,x,fi) = / 5 fc (r',x,^)exp<^ - / ct(t",x,{j,) \, (11) 

Jo A 4 I Jt' a* J 

f 7 ^ dr' - f f T ' „ dr" 1 
I k (r,x,fj,) = / S k {r ,x,-n)expi- <r(r , x, -/x) >, (12) 

J t U \ Jt M 



where a(r,x,fi) = <r cs (x) + a 77 (r, x, At). Using equation @ the source function can be written as: 

Sfc+i = R/fc • (13) 

This procedure gives the dependence of the Stokes parameters on frequency, angle and optical 
depth. Iterative methods where the calculations of the spectral structure and of the angular 
polarization structure of the radiation field were separated (Sunyaev and Titarchuk 1985; Phillips 
and Meszaros 1986) fail to obtain the frequency dependence of the Stokes vectors for a given 
scattering order. Substituting r = tt into equation (11), and r = into equation (12), we obtain 
the emergent Stokes vectors. 
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3.2. Radiative Transfer in Cylinders 

In order to treat radiative transfer in cylindrical geometry we divide the cylinder into 
horizontal spatial layers and average the computed radiation field over each layer (over the radius 
and the azimuthal directions) leaving only the dependence on the zenith angle. To simplify 
the calculations we assume the soft (reprocessed and internally produced in the cold disk) and 
reflected radiation to enter uniformly at the base of the cylinder. Thus, we effectively convert the 
2D-problem into a ID-problem. 

The boundary conditions are the same as in the slab case, but the radiation incident on the 
base of the cylinder is now: 

I iQ (x, fi) = gl re {i(x, y) + c hh B x (T hh ) + c disk B x (T disk ). (14) 

The parameter g is the fraction of the reprocessed and reflected radiation from the cold disk that 
enters the active region. In the case of slab geometry, g = 1. For cylinders atop a cold disk g ~ 0.6 
if the vertical rr equals the radial tr = Rn e aT, while g ~ 0.45 for rp = 2tr. The parameter g is 
smaller for active regions detached from the cold disk. The normalization constants, Cbb and Cdisk? 
are given by: 

/■oo roo 

c hh TT dxB x (T hh ) = gl Icpi /ir; c disk Tr dxB x (T disk ) = l disk /ir. (15) 
Jo Jo 

The expressions connecting the radiation field inside the cylinder with the source function are 
analogous to equations ( |TT| ) and (|i~2|): 

l£(r,x,y) = / S k (r',x,y)ex.p\- a(j",x,y) \C\{t,t',h), (16) 

JO {!> I Jt' A* J 

/tt dr' ~ f A 7 "' dr" 1 
Sk{T,x,-fj,)exp i -J cr(r",x,-/x) \ C 1 _{t,t' , y) , (17) 



where the correction factors, Cj_, reduce the contributions from the source function at r' to the 
radiation field at r as compared to the slab case (note that fi > 0): 

[0, if t > 1, 

c ±^ T >ri = iL rccost - tVT ^), if a<i, (18) 



t = ^I-^t-A/^tk). (19) 
The equation for the source function (O) remains unchanged. 



The emerging (polarized) flux consists of two parts: first, the radiation emerging through 
the top and the bottom of the cylinder, and, second, the radiation emerging through the vertical 
surface. The first part can trivially be found from equation ( |16| ) by multiplying the emerging 
intensity with fnr (tt appears because of the definition of the compactness). The second part is 
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given by: 

F^ de > + (x,») = - Pdr f — S k (r',x^)eK P {- [ a( T \x^) — ) Cf(r,r' , fi) , (20) 

TR Jo Jo V I Jr> V J 

Ff*>-(x,n) = — Hdr H— Sktf^-riexpl- [ T a(T",x-») — \ C*(t,t' , 

TR Jo Jt V { Jt V J 

where 

C|(r,r',/i) = 2^1 -/iVl-i 2 , if * < 1, (22) 
and equal to zero otherwise, f is given by equation (|l9l). 



3.3. Radiative Transfer in Hemispheres 

In the case of hemisphere geometry, we average over the horizontal layers just as we did for 



cylinder geometry. The incident radiation is given by equations ( |i4| ) and (15), where the parameter 
g ~ 0.7 for hemispheres atop the cold disk. The expressions ( |T6| ) and ( p"7| ) connect the radiation 
field inside the hemisphere with the source function with correction factors Cj_ being given by: 

0, if t > r + r', 
C\{tS,h) = { 1, if t<r'-r, 

C 1 , if r' — r < t < r + r', 

0, if t>r + r', 



(23) 



(r'/r) 2 , if t < r — r', 
C 1 , if r-r'<t<r + r', 



(24) 



where 



and 



^l- /j, 2 \t-t'\/ij,, r = sji 



cos 



r 2 d>* + r' 2 cos 1 



r' 2 + 1 2 " 



J2 



r 2 + t 2 



2r't 



rt sin ( 



2rt 



(25) 



The emerging flux through the base of the hemisphere is (j,ttI^(t = 0, x, fj). The flux through 
the curved hemisphere surface is given by equations (|20| ) and (21), where C± are (note that 
tt = tr and \x > 0): 

0, if t>r + r', 

2itht/t t , if t<r' - r, (26) 
C+, if r'-r<t<r + r', 

0, if t > r + r', or t <r — r', 
C-, if r — r'<t<r + r', 



C*(t,t', M ) 



(27) 
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c± 



±h4>±t/t t + yl - H 2 \l 1 - (r/r T ) 2 sin </3 ± 

0*, if H > V 1 ~ ( t / t t) 2 , 

min{ 7r/2,0*}, if A 4 < V 1 = ( t / t t) 2 , 



(28) 



3.4. Isotropic Source Function Approximation 



In some applications where high accuracy is not needed and we are not interested in the 
polarization of the radiation, we can substantially reduce the computing time assuming that 
the source function, Sk(r,x,fj,), is isotropic and homogeneous for scattering orders k > 2. The 
accuracy of this approximation is discussed in § 5.2. The approximation works because for the 
optically thin coronae photons scattered more than a few times are almost isotropic and are 
distributed almost homogeneously throughout the medium. In this case, the iteration procedure 
starting from the second scattering order can be written as follows: 

dxi 



Sk+i(x) 
Jk+\{x) 



x2 \ -^rR(x,x 1 )J k (xi), 
Jo xf 

S k+1 (x)Pj(x), k>l, 



where 



1 rr T i /-i 
J k (x) = — dr- / dfiI k (T,x,n) 

TT JO 2 J-l 

is the intensity averaged over optical depth and angles, and 

R(x,xi)= [ dfi [ d/ii [Rii(x,fj,;xi,fj,i) + R u (x, fi;xi, ~m) 
jo jo 

is the angle averaged Compton redistribution function. The quantity, Pj(x), can be obtained from 
equations (|l6l) and (|l7|): 

dr" 



(29) 
(30) 

(31) 
(32) 



l rr , l r 1 du 

dr- / — 

2i n 



TT JO 



dr' exp 



(t(t",x,h) \ C\ (r,r',/i) 



+ 



dr' exp 



dr' 



cr(r",x,-/i) \ CL(t,t',h) 



(33) 



In the case of slab geometry, where C\ = 1, and at photon energies, x < 1, where pair production 
is not important and hence <r(r, x,fi) = a cs (x), this integral can be computed analytically: 



Pj(x) 



r T - 



1 



1 (\ 



r :r 



E 3 (t x ] 



(34) 



where t x = rT<T cs (^) is the frequency dependent optical depth and £"3 is the exponential integral 
of the third order. For t x <C 1, we have 



Pt(x) 



1 



ln t x + - - 7s 



(35) 
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where = 0.577216... is the Euler's constant. 

The emergent fluxes through the top and bottom of the cylinder, and through the base of 
the hemisphere, and from the slab, can be found from equations ( |l6| ) and ( |l7| ) where the source 
functions, Sf., can now be taken out from the integrals. Corresponding //-dependent multiplicative 
factors similar to Pj(x) can be computed before the iteration procedure. Similarly, the emergent 
fluxes through the side of the cylinder, and through the curved surface of the hemisphere can be 



found using equations (|20|) and (21). 



4. The Balance Equations 

4.1. The Energy Balance 

It is common for problems where pair production is important that the luminosities appear 
in the form of dimensionless compactnesses. We define the compactnesses in the following way: 

1. For an active region (a coronal cylinder or a coronal hemisphere) 

A dissipation compactness, /diss = {LdissH / R 2 ){c t t /m e c 3 ), characterizes the dissipation 
with L(ji ss being the power providing uniform heating of the active region, and R being the 
radius of cylinder, and H its height. For hemispheres we have H = R. The soft compactness 
l s = (L S H I R 2 )(aT /m e c 3 ) characterizes the soft (reprocessed plus internally dissipated) luminosity 
from the cold disk entering the active region; l c = (L C H j R 2 )(a^ /m e c 3 ) is the coronal compactness 
corresponding to the total luminosity of Compton scattered radiation and radiation emitted in the 
corona. 

2. For a plane-parallel slab corona 

A local dissipation compactness is defined as /diss = (L diss / H)ip r v /m e c 3 ) with Ldi ss being the 
power providing uniform heating of a cubic volume of size H in the slab. Similar definitions hold 
for l s and l c . 

The soft compactness, l s , consist of two parts, /disk and gl re pr, where /disk is the compactness 
of the power that is internally dissipated in the cold disk and that enters the corona, and /repr is 
the compactness of the power reprocessed by the cold disk. The parameter g is the fraction of the 
radiation reprocessed and reflected from the cold disk which enters the active region. Introducing 
the parameter d = ZdiskAdiss; we can write the energy balance equation for the cold phase as 

^s = 5^rcpr ~\~ t^diss- (36) 

If all power dissipates in the corona then d = 0. The total coronal compactness, l c , is the sum of 
the compactness dissipated in the corona, /diss; and the part of the soft and reflected compactnesses 
which is scattered in corona: 

lc = ^diss + Psc(/s + gkcs)- (37) 
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This equation represents the energy balance of the hot phase (i.e. the corona). Here p sc is the 
probability of scattering in the corona for disk photons entering the base of the active region (p sc 
is a geometry dependent function of tt); and Z re fl is the compactness reflected from the cold disk. 
Introducing the integrated disk albedo, a, which is the fraction of the incident luminosity reflected 
by the cold disk, and the anisotropy parameter, rj, which is the fraction of all coronal radiation 
(Comptonized, annihilation, bremsstrahlung, and double Compton radiation) that is incident on 
the cold disk, we can write: 

kcfi = arjl c , l iepi = (1 - a)rjl c . (38) 
The equations (|3^)-(|38|) can easily be solved for the ratios Z c /Zdiss an d 4Adiss : 

k _ l+Pscd , . 

'diss i - Pscgv 
Is _ gyp- - a) + d(l - pscgya) 
^diss i - p sc gv 

Defining the amplification factor, A = l c /l s , we obtain 



(40) 



A = 1+Pscd . (41) 

577(1 - a) + d(l - Pscgw) 

If all power dissipates in the corona we have A = 1/577(1 — a). 

We now show how to compute the parameters entering the energy balance equations from 
the solution of radiative transfer. Let us define the partial flux emergent from the corona after k 
scatterings as (note that \i > 0): 

Fk&'V) = /^!~( r = t t,^,/j), 

F fc -(x,/i) = ^(t = 0,x,»), (42) 
for slab geometry. For cylinder geometry, the expressions are the following: 

F+(x, n) = F s k ldc ' + (x, n) + 777T/+(r = T T , X, fi), 

Fk ( X >V) = F fe slde {x,n) + fiTrl^ (t = 0,x,h). (43) 

Analogous expressions hold for hemisphere geometry, but the upward flux is just 
F£(x,ii) = F^ lde,+ (x , fi) . Let us also define F^ n,+ as the emergent flux of unscattered soft 
and reflected radiation entering the corona (note that F™'~ = 0). The emergent flux of 
unscattered radiation emitted in the corona is F . The total emergent coronal flux is then given 

by Pf = F± + Y1 F k> and the total emergent flux is F ± = F^ + F±. The total reflected 

k=l 

radiation can be found by convolving of Green's matrix for Compton reflection with the total flux 
incident on the cold disk: 

roc rl A 

I re fl(x,/j)=/ dxi / G(x,/j;xi,/ii)F~(xi,/Ji). (44) 

Jx JO Ml 



-13- 



The integrated disk albedo, a, is the ratio of the reflected flux to the total flux incident on the 
cold disk: 



Io° dx Io F 0^)°^ 
The anisotropy parameter, 77, is given by 



(45) 



= S™dxfoF c (x,ij,)dn 

Io° dx Jo [ F J~( X > t*) + p c{x, /x)J d/i 
and the scattering probability, p sc , for slab geometry is given by 

Psc = 1 - Cdx/p^^^ (4?) 
/o°° da; /o Im{x,n)fjd/Jl 

with similar expressions for hemisphere and cylinder geometry, but a factor tt should then be 
introduced in the denominator. In the calculations, d is specified and g is determined in advance, 
while a, rj, and p sc are calculated from the radiative transfer results using expressions fl45|)-(p7|) , 
and the amplification factor, A, is given by equation (f4l|). 

The sum of Zd; ss and l^isk can be written as the sum of the total upward emergent flux and the 
total downward emergent flux that does not reenter the corona: 

/*oo rl 

^ss + ^disk = 2vr / dx [F+ (x, n) + (1 - g)F~ (x, /i)] d/i. (48) 
Jo Jo 

The actual value of Zdiss is determined by pair balance, but does not influence the energy balance. 



4.2. The Pair Balance 



For the range of temperatures of interest, < 2, particle-particle and particle-photon pair 
production is negligible compared to photon-photon pair production. The pair annihilation rate, 
n ann cm -3 , is uniform throughout the corona, while the pair production rate, n 77 (r) cm" 3 , 
depends on the radiation field inside the medium being largest at the center where the photon 
density is largest and smallest at the boundaries. In pair balance, the pair annihilation rate, ?^ann: 
is equal to the volume-averaged pair production rate: 



n 



77' 



where 



1 

T T 



T T 



77 



V)dr, 



for slabs and cylinders, 



77 



(49) 



(50) 



1 - (r/rr) 



n 



77 



V)dr, for hemispheres. 



3 

2tt jo 

Here the extra factor in the integrand for hemispheres accounts for the decreasing volume of 
horizontal layers with height. The expressions for 'Harm 

and n 77 (r) are given in Appendices |A.5| 

and A. 6. 
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4.3. The Iteration Procedure 

For a given coronal geometry (slab, cylinder, or hemisphere atop of the cold disk) only two 
parameters, /diss (or, alternatively, O) and Tbb, uniquely specify the simulations if all power is 
dissipated in corona (d = 0). For d ^ 0, two more parameters, d and Tdj s k, should be specified. 
The equilibrium state satisfies the energy and pair balance equations coupled with the radiative 
transfer. To find the solution we make use of an iterative procedure. For a given 0, we compute 
the Compton redistribution matrix and the cross section, the coronal emissivities (annihilation and 
bremsstrahlung), and guess the initial values for the Zdiss an d tt- The reflected spectrum and the 
pair production absorption coefficient are set to be zero. Initial values for the parameters in the 
energy balance equations are a = 0, rj = 1/2, p sc = 0. We also compute the amplification factor A, 
and the soft compactnesses, l Tepi and Zdisk- We then normalize the incident black body radiation 
using equations (P) or (|T5|) and thus obtain the incident spectrum, Ii n (x,[i), from equations (01) or 



(14). Solving the radiative transfer by expansion in scattering orders we find the radiation field 
inside the medium as well as the emergent fluxes. We then compute the rate of pair production 
and the absorption coefficient, the double Compton emissivity, the reflected spectrum, and the 
albedo, a, the parameters 77, and p sc , and the amplification factor, A. 

By comparing ra 77 with n ann we calculate a new imposed dissipation compactness 



^dfss = ^dlss V Comparing the calculated amplification factor A new with A old we choose 

the new optical depth, tt, to be smaller than the old tt if A new > A old , and larger if A new < A old . 
The change, Att, decreases by a factor of two for on each iteration when the sign A-pp changes. 
After that we start the next iteration by again solving the radiative transfer. 

The number of iterations needed to achieve an accuracy better than 1% in all equations is 
about 10. On a Sparc 20 a typical simulation takes about 5 minutes for 6 angular points, 7 spatial 
zones, and 80 frequency points. The isotropic source function approximation (see § |3.4| ) reduces 
the computing time for solving the radiative transfer problem by an order of magnitude. 



5. Comparison with Other Codes 

5.1. Comparison with Non-Linear Monte-Carlo Results 

We compare our calculations based on the iterative scattering method (ISM) with the 
corresponding results using the Non-Linear Monte-Carlo (NLMC) code by Stern (see Stern et 
al. 1995a). We made 3 test runs each for slabs and for hemispheres atop of the cold disk. We 
assume that all power dissipates in the corona (d = 0). The parameter g is calculated in the 
iteration procedure assuming that only the radiation reprocessed below the base of the active 
region actually reenters the active region. The results are given in Table 1 and are shown in 
Figure g. We find that for a given O, the optical depth, tt, is almost the same for both codes 
with the largest difference being about 5-8 per cent at small tt- The difference in the derived 
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compactnesses, Zdiss> is less than 20 per cent for slab geometry. In the case of hemispheres, /diss 
differs by 20 per cent at large tt, and is a factor two smaller at small tt due to our approximate 
treatment of the radiative transfer in hemispheres. The procedure of averaging the radiation field 
over the horizontal layers artificially increases the photon density in the active region causing 
pair balance to be reached at smaller compactnesses. For large tt, the difference is smaller due 
to smaller boundary effects. The differences can also be due to our assumed homogeneity of the 
corona, rather than using a number of zones as in the NLMC method. As shown in Figure @, both 
codes give quite similar spectral shapes for the emerging radiation for both types of geometries. 
The differences in Z^iss can be considered as small (at least for the slab case) if we remember that 
G and tt depend rather weakly on /diss- Thus, if we fix /di ss instead of G, the differences in G and 
tt will be about 2 per cent. 



5.2. The Accuracy of some Useful Approximations 

A number of approximations can be used to decrease the time needed to compute the 
redistribution function for Compton scattering. One is the isotropic scattering approximation in 



the electron rest frame (see eq. [ A15| in Appendix |A.2| ). We found that for the mildly relativistic 
temperatures considered here the Comptonization spectra computed in this approximation is very 
accurate at small energies [x < Q) but have deficits of photons at higher energies. Solving the pair 
balance in this approximation gives Zdi ss a factor of 6 larger for large tt (small Q), and a factor of 
2 larger at small tt (large Q) (see Table 1 where this approximation is denoted ISOSCAT1). 

To improve the high energy behavior of the Comptonized spectra, we used the redistribution 



function from equation ( A15 ) and the exact value for 7* from equation ( |A8|) . This approximation 
works much better, but still produces deficit of high energy photons. We found the resulting /diss 
to be 10-30 per cent too large (see Table 1 where this approximation is denoted ISOSCAT2). 
Spectra in this approximation (dashed curves) are compared with exact results (solid curves) in 
Figure ||. At low energies (x < Q) the spectra are almost identical, but the approximate spectra 
fall more rapidly at larger energies forcing the compactness to increase in order to satisfy the pair 
balance. Due to the small contribution of photons with x ~ 1 to the total energy balance, the 
optical depth differs by less than 1 per cent from the exact calculations. We conclude that this 
approximation is useful for modeling spectra of mildly relativistic pair plasmas. 

In many works, the pair production rate is computed assuming an isotropic radiation field. 
The isotropic pair production cross section, R^°(xx\) (see Appendix A.8 ), is much easier to 



compute than the angle-dependent pair production cross section, i? 77 (xxi, /x, fj,i) (see e.g. Coppi 
& Blandford 1990). In the problem at hand, the radiation field is strongly anisotropic. We 
investigated the errors caused by making the radiation field isotropic before computing the pair 
production rate (see Table 1 where this approximation is denoted ISORAD). The resulting Zdi SS 
are systematically lower than those obtained using the exact angle-dependent i? 77 , because of 
a higher pair production rate for the isotropic case. The effect is smaller at large tt where the 
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radiation field is more isotropic. Changes in the pair production rate have no influence on the 
energy balance and consequently does not change the computed spectra and optical depths. It 
does, however, influence the value of the compactness. 

In § |3.4| we presented a method of solving the radiative transfer equation when the source 
function is assumed to be isotropic starting from the second scattering order. The overall spectral 
shapes in this approximation (dotted curves) are quite similar to the exact spectra (solid curves) 
in Figure |2[ The ISOSF approximation underestimates the flux in the "edge-on" direction for 
hemisphere geometry due to the artificial isotropization of the radiation field. The differences in 
optical depth are negligible for slab geometry, but become about 3-10 per cent for hemispheres. 
The resulting /diss are 10-50 per cent too large (see Table 1 where this approximation is denoted 
by ISOSF). 

5.3. Comparison with HM93 

In Figure || we compare our calculations with the results from Fig.4a-c in HM93 for slab 
geometry. In the low energy band they are almost identical, but at high energies (x > 0) the 
spectra of HM93 have too sharp cutoff's reflecting their use of an ad hoc exponential cutoff 
exp(— x/&). The actual cutoff energy is approximately x ~ 20, and the cutoff is not a true 
exponential, but rather reflects the thermal Compton scattering kernel and the distribution of 
the emergent photons over the scattering orders. The spectral indices in the 2 — 18 keV range, 
a2-i8) are ver y close (see Table 2) for small 0, but differs at larger 0. The reason probably lies 
in the treatment of the reflection from the cold disk. HM93 computed the reflection using the 
Monte-Carlo method assuming isotropic incident flux. For large 0, the anisotropy break (see Stern 
et al. 1995b) occurs in the 2 — 18 keV range. The flux incident on the cold disk is thus very 
anisotropic having its maximum along the normal to the disk. Making the flux isotropic by angle 
averaging artificially increases the flux along the plane of the disk, which has a larger probability 
for reflection (Matt 1993; Poutanen et al. 1996). The contribution of the reflection component to 
the total flux increases making the 2-18 keV spectra flatter. This explains the difference in «2-i8 
at large 0. 

The 0-tt relation obtained by our code agrees with calculations of HM93 to within a percent 
or two for small Tp ~ 0.01. The difference increases with increasing tt and at our largest tt = 0.37 
the of HM93 is about 10 per cent too large due to our spectral differences and the corresponding 
influence on the energy balance. The anisotropy factor, r/, (see Fig. 2a in HM93) agrees very well 
up to tt = 0.1. Above that the rj of HM93 is slightly too small, becoming 0.02 smaller at tt =0.37 
most likely due to our being 10 per cent smaller. Our albedo, a, is smaller (0.13 instead of 0.16) 
at small tt (see Fig. 2b in HM93). Above tt = 0.05, our a is 0.01 larger. The differences are 
likely due to our use of a fully relativistic and anisotropic treatment of the reflection. 

The comparison of compactnesses is not so easy. First of all, HM93 give the values 
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for the Compton compactness l c which is related to the total dissipation compactness as 
^diss = (1 — Psc?])lc, where p sc is the probability of scattering in corona (see § |4.1j ). Second, they 
define the compactness to be factor of tt larger than our definition. Thus, in Table 2 we give the 
corrected /diss = (1 — PscV)^/^, instead of the original values of l c from Table 1 in HM93. For 
small tt, our compactnesses are a factor 5 smaller, but for tt = 0.2 our compactness is larger 
than corresponding compactnesses of HM93. The approximate estimates by HM93 of the pair 
producing photon density inside the slab (see Appendix B in HM93) and the prescription for the 
Comptonized spectra (Zdziarski 1985) used in the pair balance calculations are responsible for the 
remaining differences. 



5.4. Comparison of Comptonized Spectra with Analytical Formulae 

We compare the angle averaged Comptonized spectra computed using the ISM code with 
the analytical formulae for thermal Comptonization from the papers by Titarchuk (1994, his eqs. 
[35] and [44]) and Hua & Titarchuk (1995, their eqs. [9] and [10]). We consider monochromatic 
incident photons with Iivq = 8 eV on the lower boundary of the slab. The only process which 
is taken into account is Comptonization. No pair or energy balance is imposed. Calculations 
for three different optical depths and temperatures are presented in Figure All three cases 
correspond to regime 2 in Hua & Titarchuk (1995). We find that the Hua & Titarchuk formulae 
rather well represent the general spectral behavior, but give systematically fewer photons in the 
high energy tail. The Titarchuk formulae, however, give a very good description of the spectra for 
relatively small temperatures (0 < 0.2), while they produce too many photons in the Wien bump 
for larger temperatures. In all these cases we should, however, remember that spectra below the 
anisotropy break as well as the high energy tail have quite different behavior at different viewing 
angles. The analytical formulae do not provide this angular dependence of the spectrum, and are 
therefore quite limited in practice. 



5.5. Polarization Properties 

In this section, we compare our calculations of the degree of polarization of the radiation 
emerging from the Compton scattering slab-corona with some earlier results by others. We 
assume here that the cold disk emits semi-isotropic unpolarized radiation. We define the degree 
of polarization as p = (Q/I) 100%. The polarization is positive when the electric vector is 
predominantly parallel to the normal to the slab. The behavior of the total polarization is affected 
by both the Compton scattering radiation from the hot corona, and by reflected radiation from 
the cold disk. 

Sunyaev & Titarchuk (1985) calculated the polarization of the Comptonized radiation from a 
slab. The angular and polarization structure of the radiation field were obtained with an iteration 
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procedure based on an expansion in scattering orders using the Rayleigh matrix. The frequency 
dependence of the intensity was obtained solving the Kompaneets equation (see Sunyaev & 
Titarchuk 1980), which is a diffusion equation in frequency space. The distribution of photons 
over escape time and thus the relation between the scattering order and frequency were obtained 
by solving the diffusion equation in physical space. This approach has a few shortcomings. First, 
the diffusion approximation in optical depth (i.e. physical space) is not valid for the optically thin 
coronae considered here. Second, for the case of a hot electron gas, the frequency redistribution 
cannot be considered as diffusion due to the large frequency shift in each scattering. And last, 
for high electron temperatures (G <; 0.1) and/or large photon energy (x <L 0.1), the polarization 
properties of the exact Compton redistribution matrix are quite different from those of the 
Rayleigh matrix (Poutanen & Vilhu 1993). Thus, even for quite small electron temperature, 
G = 0.11, and relatively large optical depth, tt = 0.5, the maximum polarization of the hard 
radiation computed by this method is p « 50% (see Fig. 8 in Sunyaev & Titarchuk 1985, note also 
that their To = tt/2), compared to p « 25% in our calculations (see our Fig. ||). The differences 
become much larger for smaller tt and/or larger Q. 

Following the same idea of separating the polarization structure from the frequency 
redistribution, Haardt & Matt (1993) computed the degree of polarization from an optically thin 
hot slab-corona applying the method of Haardt (1993) to obtain the spectra. Using an iterative 
scattering scheme similar to ours they avoided the diffusion approximation both in frequency 
and optical depth space. They, however, still used the Rayleigh matrix. To compare the results, 
we have chosen the same parameters as in Figure 2 in Haardt & Matt (1993). The blackbody 
temperature is taken to be, Tbb = 10 eV, the optical depth is tt = 0.5 and 0.05. We show the 
results of our computations in Figures || and ^. We present the results for two viewing angles, 
H = 0.11 and [i = 0.50. The degree of polarization is zero in the direction normal to the slab due 
to symmetry. It is clearly seen that the degree of polarization is much too large in Haardt & 
Matt (1993). They found the maximum polarization to be p ~ 45% for tt = 0.5 and p ~ 33% for 
tt = 0.05, while we obtained p « 25% and p « 5%, respectively (thick solid curves in Figs. || and 
!)• 

Below, we discuss the polarization properties both of the Comptonized radiation from the hot 
corona and of the radiation reflected from the cold slab. In order to better see the contribution 
to the overall polarization from different scattering orders we also show the spectra for individual 
scattering orders as well as the total spectrum in the upper panels of Figures || and |6[ 

First we consider the polarization caused by scattering in the hot corona. The degree of 
polarization increases with the number of scatterings reaching its asymptotic value after few 
scatterings. The asymptotic value depends on optical depth, temperature and zenith angle. 
At energies close to m e c 2 , the Klein-Nishina corrections start to be important decreasing the 
polarization. As the electron temperature increases the polarization decreases (Poutanen & Vilhu 
1993). At a given frequency the contribution from the higher scattering orders (with larger 
polarization) becomes smaller also decreasing the polarization. Due to these two reasons the 
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polarization in Figure || (tt = 0.5) is larger than in Figure |6| (tt = 0.05). The polarization 
of a given scattering order first decreases toward higher energies but then increases due to the 
contribution from the scattered reflected component. The decrease is caused by the fact that 
the largest energy shift is obtained in backward scatterings which do not produce additional 
polarization. 

The reflection component contributes to the overall spectra at energies x ~ 0.01 — 1. When 
the flux incident on the cold disk is nearly isotropic, which is the case for optically thin corona, 
the polarization of the reflected component is positive at x ^ 1. It is maximal in edge-on 
directions, and is zero in the normal direction (see Poutanen et al. 1996 for a discussion of the 
polarization properties of Compton-reflected radiation from the cold slab). The polarization in 
the direction close to the normal decreases at higher energies, changing sign at x ~ 0.6. The flux 
reflected close to the normal direction cuts off at energies above x ~ 1. At directions closer to the 
plane of the slab, the cutoff is slower, and is determined by the cutoff of the incident spectrum. 
The polarization has a sharp feature at 6.4 keV due to the contribution from the unpolarized 
fluorescent iron line. 

For tt = 0.5 (Fig. ||) the polarization of the reflected radiation (and that of the reflected 
component scattered once or twice in the corona) is smaller than the polarization of the component 
produced by Compton scatterings in hot corona. This causes the polarization to decrease at 
x > 0.1. On the other hand for tt = 0.05 (Fig. ||) the reflected component has significantly larger 
polarization than the scattered component resulting in a smoothly increasing polarization in the 
energy interval from x ~ 0.01 up to x ~ 0.1. The decrease of the polarization of the reflected 
component at higher energies and the change of sign at x ~ 0.6 for [i = 0.5 cause the drop in the 
total polarization at x ~ 1. 

Future observations of X-ray polarization by Spectrum-X-'j satellite (Kaaret et al. 1992) can 
be a powerful tool for determining the physical conditions and, probably, the geometry of the 
X-ray emitting region in AGN and X-ray binaries. We note here that the degree of polarization 
for hemisphere and cylinder geometries is smaller than for the slab case. We can conclude that 
if small polarization in the X-rays will be observed this could argue for Comptonization models 
where the temperature of the electron gas is large and/or for models where the geometry of the 
corona is not slab-like. 



6. Summary and Conclusions 

We have described a versatile code based on the iterative scattering method (ISM) to 
accurately solve the radiative transfer and Comptonization in a two-phase disk-corona models for 
active galactic nuclei and X-ray binaries. 



The code has several attractive features some of them being unique to this code: 
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1) The radiative transfer is fully angle-dependent, and one can easily determine the outgoing 
spectrum in any direction to high accuracy at any photon energy. The outgoing spectrum in a 
given direction may differ greatly from the angle-averaged spectrum. 

2) The radiative transfer is valid for both nonrelativistic and relativistic temperatures. 

3) The radiative transfer is polarized. 

4) Most important radiation processes in hot thermal plasmas including Compton scattering, 
photon-photon pair production, pair annihilation, bremsstrahlung, and double Compton scattering 
are taken into account. The latter two were not important for the parameters of the test cases we 
considered in this paper. 

5) The corona can be in energy and/or pair balance. 

6) The corona can either be a pure pair corona, or one can include a background plasma. 

7) The reflection by the cold disk is exactly treated using a reflection matrix that, in 
particular, accounts for the full angular dependence of the incident radiation. 

8) The ISM code has been extensively tested against a Non-Linear Monte Carlo (NLMC) 
code (Stern et al. 1995b) finding very good agreement. 

9) The ISM code is an order of magnitude faster than the NLMC code. The ISM code also 
allows for an easier determination of the spectral fall-off at photon energies above kT e . The ISM 
code also gives more accurate emerging spectra at a given viewing angle as compared to the 
NLMC code where one must average over a range of viewing angles in order to improve the photon 
statistics. 

10) Various approximations for the radiative transfer/Comptonization can be used in order to 
improve computing efficiency Quite accurate results can be obtained if one assumes the Compton 
scattering source function to be isotropic. The gain in computing efficiency is then an order of 
magnitude. A typical run on a Sparc 20 then takes about 40 s for three angular, 80 frequency, and 
7 spatial gridpoints. 

There are, however, some limitations: 

1) The iterative scattering method converges only for small optical depths. For small 
temperatures (0 < 0.1) and large optical depths (r > 1), the round off errors become large due 
to the large number of scatterings and the accuracy of the results decreases. The maximum 
allowed tt depends on the temperature and is approximately 1 for 6 = 0.1, and 1.5 for = 0.5 
in the slab case and 2 — 3 in the case of hemisphere geometry. This is not much of a limitation 
at temperatures above about 100 keV, as the optical depths needed to explain observed X-ray 
spectral indices in AGN are necessarily less than unity. 

2) The ISM code is one-dimensional. The ISM code can, however, be applied to quasi-lD 
radiative transfer in two-dimensional active regions with cylinder or hemisphere geometry atop or 
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elevated above the cold disk. The NLMC code can, in principle, treat arbitrary geometries. 

3) The corona is uniform in temperature and density. The Nonlinear Monte Carlo code has 
been used (Stern et al. 1995b) to show that the differences in temperature and density across the 
corona is at most a factor 2. In principle, one could divide the corona into a few zones, and solve 
the radiative transfer equation using the appropriate Compton redistribution functions and cross 
sections corresponding to the temperature in each zone. The pair and energy balance equations 
should then be solved in each zone separately. This, however, is much more time consuming than 
considering a homogeneous corona. The resulting spectra from homogeneous and inhomogeneous 
corona are very similar. 

4) The ISM code treats steady radiative transfer. The NLMC code, on the other hand, can 
treat time dependent situations. 

There are several possible applications of the ISM code: 

1) Results of earlier work using various approximations can be checked, and the validity of 
the approximations tested. In particular, we find that the full Compton redistribution matrix, 
rather than the Rayleigh scattering matrix, must be used in order to obtain accurate polarized 
X-ray spectra. 

2) The validity of published analytical fits for angle- averaged Comptonized spectra can be 
checked. The results of the ISM code can be used to obtain analytical fits of spectra as function 
of viewing angle and for different geometries. 

3) The ISM code has great potential for modeling X-ray and 7-ray spectra from active 
galactic nuclei and X-ray binaries. The ISM code has already been used together with a NLMC 
code to interpret the statistics of observed X-ray spectral indices and compactnesses from Seyfert 
1 galaxies (Stern et al. 1995b). The anisotropy of outgoing spectra is an important ingredient in 
this interpretation, which could not have been done using angle-averaged model spectra. 

4) The ISM code is highly suitable for inclusion in spectral fitting software such as XSPEC, 
whereupon observed X-ray spectra can be modeled. Such modeling of a number of sources will 
appear in forthcoming work. 

5) The ability of the ISM code to treat polarized radiative transfer makes it a powerful 
tool for interpreting future observations of X-ray polarization from, e.g., the Spectrum-X-'j and 
INTEGRAL satellites. 

The authors thank Boris Stern for stimulating discussions. This research was supported by 
grants and a postdoctoral fellowship (J. P.) from the Swedish Natural Science Research Council. 
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A.l. The Electron Scattering Source Function and the Compton Redistribution 

Matrix 



The thermal electron scattering source function, S(t,x,/i), for an azimuth-independent 
radiation field accounting only for linear polarization can be expressed in terms of the 
azimuth-averaged Compton redistribution matrix, R(x, fi; x±, fii), as: 

S(t,x,h)=x 2 —t dniR(x,fj,;xi,m)I(r,xi,fj,i) , 
Jo x 1 J-l 



(Al) 



Here, R is the azimuth-averaged product of two rotational matrices, L, and the thermal Compton 
redistribution matrix, C(x,xi,cos6) (Poutanen & Vilhu 1993; the hat identifies R, L, and C as 
matrices, the tilde identifies S and I as vectors): 



r 2 * 

R(x,fi;x 1 ,[i 1 )= / d(pL(-x)C(x,x 1 ,cose)L(xi)- 
Jo 



(A2) 



In general, it is 4 x 4 matrix. The rotational matrices are given by the following expression (see, 
e.g., Chandrasekhar 1960): 



L(X) = 



( 1 \ 

cos(2x) sin(2x) 

-sin(2x) cos(2x) 

1 



(A3) 



Due to azimuthal symmetry and the absence of circular polarization we consider only the 2x2 
matrix in the upper left corner of the general matrix: 



R 



Ru R12 
R21 R22 



(A4) 



The elements of this matrix are: 

R u = J Cd^, 

R12 = J Ci cos 2xidip, 

R21 = J Cicos2xd<£, 

R22 = J[C+ cos 2( X -Xi) + C- cos 2(x + xi)]d¥>, 

where C± = (Cq ± C\j)/2, and the cosines are given by 

ui — a cos 6 11— iii cos 9 

cos X = . - 7= , cosxi 



sin 0y/\ — fi 2 ' 



sin6*Wl — ji{ 



(A5) 



(A6) 
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with cos 9 = im\ + \J\ — fi 2 J 1 — fj, 2 cos ip. Finally, the functions C, Ci,Cq, and C\j are four 
of the five functions forming the thermal Compton redistribution matrix, C(x, x\, cos 9). This 
matrix is given by a single integral over the electron energy distribution /(j) and the Compton 
redistribution matrix, C m (x, x±, cos 9, 7), for an isotropic monoenergetic electron gas with Lorentz 
factor 7 (Nagirner & Poutanen 1993): 



C(x, 27, cos 9) 



( c c x \ 

Ci C Q 

Cu 

y Cy J 



f^)C m (x,x 1 ,cos9, 1 )d 7 



where 



/(7)d7 



7* 



xi + Q(l + 2/q) l l 2 ]/2, g = xx t (l - cos9), 



( c m 


c m 





\ 


rjm 


rim 

C Q 














rim 













rim 
°V / 


Q 2 


= {x 


-XI 


) 2 + 2g. 



(A7) 



The function C m is the scalar redistribution function derived by Jones (1968) (see also Coppi 
& Blandford 1990). We use the following expressions to calculate the five functions that enter 



C m (x, x±, cos 9, 7) in equation (|A7|) : 



rim 
°U 

rim 

rim 
W 



rim 1 /-im 
°a b ' 

c™ + c™, 

2 u- 
Q rq 

rim I rim 



Q 



'U 

'h 



<?c a m , 







rq 



(2Q + u) - 4 



+ -+2E- 



(A8) 



where 



rim 



CI: 



(u 2 - Q 2 ){u 2 + 5v) Q 2 
u Q h it 




9 9 * 



(A9) 



and 



it = ai 
(7 - x) 2 + r, 



a? 



(x + xi)(27 + aci — x)/ (a + ai) 



aai, 



(A10) 



(7 + x\) 2 + r, 



(1 + cos0)/(l - cos( 
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The normalized relativistic Maxwellian distribution is given by 

e -7/e 

> W = 4^(1/6) - < AU) 

which gives the density of particles in the dimensionless momentum volume, 4irz 2 dz, normalized 
to unity. Here, K n is the modified Bessel function of second kind of order n, and z = \J-f 2 — 1. 
Methods for computing the integrals when averaging C m (x, x±, cos 6, 7) over a relativistic 
Maxwellian electron distribution are given in Poutanen (1994). If the electron temperature is not 
very high (6 < 1) we can use Gauss-Laguerre quadrature. 



A. 2. The Compton Redistribution Function for Isotropic Scattering in the Electron 

Rest Frame 



A very simple expression for the redistribution function, C m , can be derived if we assume that 
the characteristic photon energy in the electron rest frame is small X17 <C 1, i.e. the Thomson 
limit. The scattering in this limit can be assumed to be coherent in the electron rest frame. We 
also assume that the scattering is isotropic in the rest frame. Such simplifications give a correctly 



normalized redistribution function at energies x <C G (see eq. [A23]), whose shape slightly differs 
from the exact one. The function C m in that approximation becomes (Arutyunyan & Nikogosyan 
1980) 

4 

6 "3Q' 

being non-zero when 



(A12) 



2<?( 7 2 - 1) > (xi 



x) 



(A13) 



Integration over cos 9 between the limits defined by equation ( |A13| ) gives the angle-averaged 
redistribution function (Rybicki & Lightman 1979): 



C m (x,x 1)7 ) 



C m sinfldfl 



3xx\ 



■ 7 I 
X + X\ \x 

z 



Xl 



(A14) 



where x/x% G [(7 — z) 2 , (7 + z) 2 ]. For Maxwellian electrons the redistribution function is given by 



C(x, Xl, cos l 



-7*/e 



87rQK 2 (l/G)' 



(A15) 



where 7* = Q/y/2q. Equation ( |A15 ) gives very good approximation to the exact redistribution 
function for mildly relativistic temperatures and x <C G. Notice also that xC(x,x\, cos 6) is a 
function of the ratio xjx\. In Table 1 this approximation is called ISOSCAT1. Even better 
agreement with the exact redistribution function is obtained if 7* from equation (A8) is used 
in equation ( |A14 ). We call this approximation ISOSCAT2 in Table 1. The accuracy of these 
approximations is discussed in § p\2 
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A.3. The Thermally Averaged Compton Scattering Cross Section 



The Compton scattering cross section averaged over a relativistic Maxwellian electron 
distribution can be written as an integral over the electron energy: 

1 + 2x( 7 + z) 



3<7t 



16x 2 6^ 2 (l/6) 
2 



-f 

9)7i 



-7/e 



9 2 
a?7 + - + -7 ) In 

2 x 



+ z \ x — 



ln(l + 4x7 + 4 ^ ) + 



1 + 2x(7 - z) 
4x 2 z(7 + x) r^ii+z) 
1 + 4x7 + 4 ^ 2 Jx(y-z) 



2xz 



- 2 / ln(l + 2^ 

Jx(-/—z) 



d£ 



d7(Al6) 



Making the substitutions 7 = 1 + 6 exp(— 2t) on the interval [1, 1 + 6], and 7 = 1 + 6(1 + t) on 
the interval [1 + 6,00), and applying 10-points Gauss-Laguerre quadrature formula we achieve an 
accuracy better than 0.02 per cent. 

In limiting cases, the thermal cross section can be computed using simple expressions (Gould 
1982; Svensson 1982; Nagirner & Poutanen 1994): 



0" cs (x) 
(T cs {x) 



3ut 
8x2" 

3d"T 

16x6 

0t 



4+ (x-2- ^ ln(l + 2x) + 



2x 2 (l + x) 



- -7s + ln4x6j , 



, V(-&)v„ +2 (i/e), 

^2(1/6) n=0 



where 



n + 2 + 



+ 



(l + 2x) 2 
6 > 1, x6 > 1, 

x6 < 1, 



16 



e< 1, 



(A17) 
(A18) 
(A19) 

(A20) 



n+l n + 2 ra + 3 

These simple approximations can be used to check the correctness of the thermal cross section 
routine and to estimate the numerical accuracy. 



A. 4. Symmetry Properties and a Normalization Condition 

The azimuth and thermally averaged functions Rij have symmetry properties which can be 
exploited to simplify the radiative transfer equation, to reduce the time needed to calculate all 
elements of the redistribution matrix, and to check the accuracy of the calculations: 

1. Frequency symmetry 

i^(x,/i;xi,^i)e _Xl/e = R ij (x 1 ,iJ,;x,ii 1 )e- x/e , i,j = 1,2, (A21) 

which follows from microscopic detailed balance between states x and xi when the photons and 
the electrons have a Wien and a Maxwellian distribution, respectively (Pomraning 1973; Meszaros 
& Bussard 1986). The exponential factors represent the Wien distribution, while the photon phase 
space factors have been absorbed in the definition of Rij. 
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2. Angular symmetry 

Rij(x,fj,;xi,fj,i) = Rji(x,fj,i;xi,/j,), 

Rij(x,fj,;xi,fj,i) = R ij (x,-fi;x 1 ,-fi 1 ), i,j = 1,2. (A22) 

These angular symmetries follow directly from the fact that the scattering process depends on the 
scattering angle between the in and outgoing photons, and not on their angle cosines, fx and 
separately. 

The thermal Compton scattering cross section, cr cs (x), and the scalar redistribution function, 
Ru(x±, fi±; x, fi), i.e. element 11 of the thermal redistribution matrix, are related through a 
normalization condition (see e.g. Pomraning 1973; Nagirner & Poutanen 1994): 

a (x) 1 f°° f 1 f 1 

— = -/ xidxi dfi dfii[Rn(x 1 ,fj, 1 ;x,fi)+Rn(x 1 ,fii;x,-fi)]. (A23) 

ct x Jo Jo Jo 

Analogous integrations of R12 and R21 gives zero on the left hand side due to the requirement 
that the polarization is zero (Q = 0) for an isotropic radiation field. These relations can be used 
to check the accuracy of the calculation of the redistribution matrix, and to estimate the quality 
of the frequency and the angular discretization. 



A. 5. The Photon-Photon Pair Production Rate 

For azimuth-independent (i.e. axisymmetric) photon distributions, the rate of photon-photon 
pair production, n 77 cm -3 s _1 , neglecting polarization is given by integrals over dimensionless 
photon energy, dx, and solid angle, 2ndfi, as follows: 

. . 1 2vr f°° dx f 1 T . . . . . 

n nK T ) = o 2/ — / I{r,x,n)a^(T,x,n)dii. (A24) 

2 m c c 2 Jo x J -1 

Here, 7(r, x, /x)dx/(cm e c 2 x) is the number density of photons of energy x in the interval dx 
per steradian traveling in the direction /j,, and a 77 cm" 1 is the absorption coefficient due to 
photon-photon pair production. The factor 1/2 is due to both interacting species being photons. 
The absorption coefficient is given by another integral over the target photon energy, dxi, and 
solid angle, d(pdfj,±, 

r 2 f°° dx\ f 1 

a 77 (r,x,^) = — / R-f-f(x,fjL;xi,(ii)I(T,xi,tii)diJ,i, (A25) 

m c r Jo x\ J -1 

where r e is the classical electron radius, i? 77 (x, \i\ x\, fi±) is the dimensionless, azimuth-integrated 
pair production cross section, 

/*7T 

Ryj(x, /x; xi, jJL\) = 2 / s 77 (cj)(1 — cos 0)d(p, (A26) 
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the (1 — cos 9) factor is discussed in, e.g., Weaver (1976), and s 77 (u;) is the dimensionless 
photon-photon pair production cross section (Jauch & Rohrlich 1976) 



2 + - - -k) cosh" 1 ^- + 



(A27) 



Here, (1 — cos#)/2, where x cm is the photon energy in the center-of-momentum 



frame, and 9 is the interaction angle related to other cosines as cos 9 = fi/ii + — /i 2 y 1 — ^\ cos V 5 - 
Using these relations, the threshold condition for pair production, x cm > 1, can be 
written as a constraint on cos 9 or on cos tp giving a minimum allowed value for ip: 



cos(/? m i n = (1 — fifii — 2/xx\)/ yV^ ~ A 4 y 1 ~~ Mi )■ The factor 2 in equation ( A26 ) comes 
from the integration range originally being </7 m i n to 2ir — </? m i n , and the integrand being an even 
function of ip around tp = ir. The axisymmetric pair production rate was previously considered 
by Stepney & Guilbert (1983) who chose x cm as integration variable instead of (p. Their rate is a 
factor 2 too large as pointed out by Kusunose (1987). 



A. 6. The Thermal Pair Annihilation Rate 

For a relativistic Maxwellian electron (and positron) distribution, the pair annihilation 
reaction rate, n an n cm -3 s , can be written as a one-parameter (0) single integral (Weaver 1976). 
Svensson (1982) made a simple fit to that integral accurate to within 2 per cent: 

„2 7T 



n ann = «-» + «-e 1 + 2e2/ln(L3 + 2 ^ er (A28) 



where r/E ~ 0.5615. 



A. 7. Pair Annihilation Emissivity 

The emissivity due to thermal pair annihilation, e ann (x, 0)dx erg cnr 3 s _1 sr _1 , in an energy 
interval dx can be written using detailed balance arguments in terms of the pair production cross 
section in the following form (Svensson 1983): 



(x, G) = n_n + r e 2 m e c 3 27re ^ 2(i/e) / ujs^e^du;, (A29) 



where s 77 (cc>) is given by equation (A27). Simple analytical fits for the one-parameter (x0) integral 
in equation ([A29| ) accurate to within 0.04 per cent are given by Svensson, Larsson, & Poutanen 
(1996). 
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A.8. Symmetry Properties of the Azimuth-integrated Pair Production Cross 

Section 

The azimuth-integrated cross section, i? 77 , obey energy and angular symmetry relations, 
which are useful in reducing the computing time: 

R 71 (x,n;x 1 ,fi 1 ) = i? 77 (xxi,/i,/ii), 

Ryy{x, /i; xi, pL\) = ii 77 (x, x±, fi), (A30) 
Ry-y(x, /i; x\, fix) = i? 77 (x, — |U; xx, — Hi). 



These symmetry relations follow directly from equation ( A26| ), the definition of u, and the 



relations for cos# and cos</? m ; n . To check the accuracy of our calculations we integrate the 
azimuth-integrated cross section over one cosine angle and average over the second in order to 
obtain the fully solid angle-integrated cross section for the isotropic case (which is a well known 
function first computed by Gould and Schreder 1967): 

R 1 *°(xx 1 ) = ^j dfi J dfM 1 H yi (x,iJ l ;xi,ni) = J dfi J d/^ [Ry 7 (x, fj,; x x , n x ) + Ryyix, fi; x u —Ml)] 

(A31) 

Here we have used the third symmetry property above. The angle- averaged function, (p{xxx) in 
Gould & Schreder (1967) is related to our R™°(xxx) through <f>{xx\) = (xx 1 ) 2 R^(xx 1 )/8tt 2 . The 
angle-averaged cross section, R(xxx) in Coppi & Blandford (1990), is related to our R^°(xxi) 
through R(xxx) = cr 2 R™(xxx) / '(4tt) . Coppi & Blandford (1990) give a useful fit for R(xxx) 



accurate to within 7 per cent for all xxx- We find our R^°(xxx) computed using equation (A31) 
to typically be accurate to within 2 per cent. 

One can show that the annihilation emissivity can be written as an integral, not over the pair 
production cross section, but over the angle-integrated pair production cross section: 

^.3g— a;/© i-oo 

earing, 6) = n„n+r 2 m e c 3 ^292^2^/9) j Q x 2 R™{xxx)e- Xl/e dxx. (A32) 

By numerically computing this integral using our computed R^°{xxx) and comparing with 
equation ( A29| ) we obtain an extra check of the consistency of our pair production and annihilation 
routines. 



A. 9. Double Compton emissivity 

The angle-averaged double Compton spectral emissivity, ejic(x,®)dx erg cm' 3 s _1 sr _1 , in 
an energy interval dx is given by the expression (see, e.g., Svensson 1984): 

t n\ 1 j_ ^ e ~ X/0 f°° -3w u f°° da DC( s ( u/x 1 + x 1 /u \ 
eDC (x,G) = (n + + n_)x^ IMj / o x x J{xx)dx x J q u—(x,u)ex V ( j do, 

(A33) 
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Here J(x\) is the mean intensity of the interacting photons. The differential cross section for the 
double Compton process is given in Svensson (1984, eqs. [A5], [A6]). In order to account for the 
high energy cutoff at x > we introduced an ad hoc exponential factor, e~ x / & . 



A. 10. Bremsstrahlung emissivity 

The emissivities due to relativistic electron-electron and positron-positron thermal 
bremsstrahlung, e±±(x, @)dx erg cm~ 3 s _1 sr -1 , in an energy interval dx are given by the 
expression: 

2 3 /2 

e±±{x, 9) = n^o-Tafm^e-^e- 1 / 2 — g±±{x, 0). (A34) 

A similar expression holds for the electron-positron emissivity e_] (x,0). We used the 

approximations by Skibo et al. (1995) for the Gaunt factors g±±(x,Q) and (x,Q). Note that 

the term of unity in equations (A7), (A9), and (A13) in Skibo et al. (1995) should be deleted. 



A. 11. Numerical Integration 



All azimuthal integrations are made using an 11-point Simpson quadrature. Furthermore, we 
apply 3-point Gaussian quadrature to calculate integrals over zenith angles for each hemisphere. 
The integration over frequencies is performed using rectangular quadrature on a logarithmic 
frequency scale, dx/x = din a;, with bin width ~ 0.1. The integrals over optical depth are 
calculated using rectangular quadrature. The number of points, N T , is dependent on the geometry 
and the optical depth. We typically used N T = 6 for slabs, and N T = 11 — 21 for active regions. 
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Table 1. Comparison of results obtained using the ISM and NLMC codes. 







SLAB 






HEMISPHERE 




Method 


e 


TT 


^diss 


e 


T T 


^diss 


ISM 


0.19 


0.29 


260 


0.24 


0.70 


560 


NLMC 


0.19 


0.29 


300 


0.24 


0.70 


700 


ISOSCAT1 


0.19 


0.29 


1830 


0.24 


0.70 


2730 


ISOSCAT2 


0.19 


0.29 


390 


0.24 


0.70 


740 


ISORAD 


0.19 


0.29 


250 


0.24 


0.70 


470 


ISOSF 


0.19 


0.30 


340 


0.24 


0.76 


640 


ISM 


0.29 


0.17 


19 


0.49 


0.29 


27 


NLMC 


0.29 


0.16 


20 


0.49 


0.28 


40 


ISOSCAT1 


0.29 


0.17 


86 


0.49 


0.29 


70 


ISOSCAT2 


0.29 


0.17 


25 


0.49 


0.29 


32 


ISORAD 


0.29 


0.17 


17 


0.49 


0.29 


20 


ISOSF 


0.29 


0.17 


23 


0.49 


0.30 


31 


ISM 


0.82 


0.036 


0.24 


1.20 


0.073 


1.5 


NLMC 


0.82 


0.033 


0.20 


1.20 


0.070 


3.0 


ISOSCAT1 


0.82 


0.037 


0.47 


1.20 


0.076 


2.4 


ISOSCAT2 


0.82 


0.036 


0.27 


1.20 


0.073 


1.6 


ISORAD 


0.82 


0.036 


0.16 


1.20 


0.073 


0.8 


ISOSF 


0.82 


0.036 


0.27 


1.20 


0.075 


1.6 



Input parameter: O = kT c /m c c 2 , the dimensionless pair temperature of the corona (volume averaged in the 
NLMC case); Output parameters: tt, the (averaged) Thomson optical depth; /diss = (£diss / H)(ctt /m e c 3 ), the local 
dissipation compactness of the coronal slab or hemisphere. In all cases, d — and Tbb = 5 eV. ISM and NLMC 
represent results using the ISM and NLMC codes, respectively. ISOSCAT1 and ISOSCAT2 represent ISM results 
using isotropic scattering in the electron rest frame (Appendix A2). ISORAD represent ISM results when the radiation 
field is made isotropic before solving the pair balance. ISOSF represent ISM results when the source functions, Sk>2, 
are assumed to be isotropic and homogeneous (see § 1.4). 
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Code 


e 


' 1 


^diss 




f>Q 1 Q 

^z— lo 




TQA/T 
lOlVl 








,. n 1 1 q 


fi — u.o 


,, n 887 

H — U.oo / 


HM93 








H = 0.15 


H = 0.55 


H = 0.95 


ISM 


0.26 


0.19 


34 


1.15 


1.11 


1.08 


HM93 


0.26 


0.20 


22 


1.23 


1.14 


1.10 


ISM 


0.67 


0.05 


0.49 


0.89 


0.77 


0.69 


HM93 


0.67 


0.05 


2 


0.87 


0.73 


0.65 


ISM 


1.75 


0.01 


0.02 


1.37 


1.16 


0.85 


HM93 


1.75 


0.01 


0.1 


1.23 


0.95 


0.53 



Table 2: Comparison of ISM results for coronal slabs with those of HM93 



Note. — Input parameter: O = kT c /m c c 2 , the dimensionless pair temperature of the coronal slab; 
Output parameters: tt, the vertical Thomson optical depth of the coronal slab; 
idiss = (I/diss/i7)(o"T/m c c 3 ), tne dissipation compactness of the coronal slab; 
CK2-18, the least square fitted 2-18 keV intensity slope for three specified cosine angles. 
In all cases, d = and T bb = 5 eV. See text for determination of ldi BB from HM93. 
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2 




log x 

Fig. 1.— (a) Emergent flux, xF x (arbitrary units), from slab pair coronae as function of 
dimensionless photon energy, x = hu/m c c 2 , for the parameters in Table 1. The spectra are averaged 
over viewing angles 0.6 < fi < 0.9. Solid and dotted curves: results using the ISM code and NLMC 
code, respectively, (b) Same as (a), but for hemisphere coronae. 
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log x 

Fig. 2. — Same as Figure [l] but for two viewing angles, \i = 0.11 and = 0.89. The face-on 
spectra (fi = 0.89) can be identified by their more prominent black body component. All curves 
represent results from the ISM code. Results using the exact scattering kernel are shown by 
solid curves. Dashed curves represent results assuming isotropic scattering in the electron rest 
frame (ISOSCAT2, see Appendix A2). Dotted curves correspond to the isotropic source function 
approximation (ISOSF). Notice that the emergent spectra for the different approximations are 
normalized to the corresponding value of /diss in Table 1. The spectral fluxes in a given direction 
differ slightly, but the spectral shapes are almost identical. 
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Fig. 3. — Emergent flux, xF x (arbitrary units), from slab pair coronae as function of dimensionless 
photon energy, x = hv/m c (? . Solid curves: results using the ISM code for [i = 0.11 and [i = 0.89. 
Dotted curves: results from Figures 4a,c in HM93 for ^ = 0.15 and fi = 0.95. 
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log x 

Fig. 4. — Angle averaged Comptonized spectra, xF x (arbitrary units), from slab coronae as function 
of dimensionless photon energy, x = hu/m c c 2 . Soft monochromatic photons with huQ = 8 eV are 
assumed to be incident on the lower boundary of the slab. Solid curves: results using the ISM code; 
dotted curves: spectra using the analytical formulae (35) and (44) from Titarchuk (1994); dashed 
curves: spectra using the analytical formulae (9) and (10) from Hua & Titarchuk (1995). No pair 
or energy balance is imposed. 
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Fig. 5. — Emergent flux, xF x , from slab coronae (upper panels) and the degree of polarization, p 
(lower panels) as function of dimensionless photon energy, x = hu/m c c 2 , for = 0.11, r T = 0.50, 
Tbb = 10 eV for two viewing angles \i = 0.11 and fi = 0.50. No pair or energy balance is 
imposed. Thick solid curves: the overall spectra and polarization. The contribution from some of 
the scattering orders is also shown. The labels show the order of scattering. The zeroth (dotted) 
component consists of the unpolarized blackbody disk radiation as well as the radiation reflected 
from the cold disk. The scattered components consist of multiple scattered blackbody radiation 
and multiple scattered reflected radiation with the latter being centered around x ~ 0.1. 
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